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Abstract 

We study the transport properties of nonautonomous chaotic dynamical systems 
over a finite time duration. We are particularly interested in those regions that remain 
coherent and relatively non-dispersive over finite periods of time, despite the chaotic 
nature of the system. We develop a novel probabilistic methodology based upon trans- 
fer operators that automatically detects maximally coherent sets. The approach is very 
simple to implement, requiring only singular vector computations of a matrix of tran- 
sitions induced by the dynamics. We illustrate our new methodology on an idealized 
stratospheric flow and in two and three dimensional analyses of European Centre for 
Medium Range Weather Forecasting (ECMWF) reanalysis data. 

1 Introduction 



Finite-time transport of time-dependent or nonautonomous chaotic dynamical 
systems has been the subject of intense study over the last decade. Existing 
techniques to analyse transport have evolved from classical geometric theory of 
invariant manifolds, where co-dimension 1 invariant manifolds are impenetrable 
transport barriers. In this work we take a very different approach, based on 
spectral information contained in a finite-time transfer (or Perron-Frobenius) 
operator. Our technique automatically identifies regions of state space that 
are maximally coherent or non-dispersive over a specific time interval, in the 
presence of an underlying chaotic system. These regions, called coherent sets, 
are robust to perturbation and are carried along by the chaotic flow with little 
transport between the coherent sets and the rest of state space. Thus, these co- 
herent sets are ordered skeletons of the time-dependent dynamics around which 
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the chaotic dynamics occurs relatively independently over the finite time con- 
sidered. We develop the theory behind an optimization problem to determine 
these coherent sets and describe in detail a numerical implementation. Numer- 
ical results are given for a model system and real-world reanalysed data. 

Transport and mixing properties of dynamical systems have received considerable interest 
over the last two decades; see e.g. |l|-l4j for discussions of transport phenomena. A variety 
of dynamical systems techniques have been introduced to explain transport mechanisms, to 
detect barriers to transport, and to quantify transport rates. These techniques typically fall 
into two classes: geometric methods, which exploit invariant manifolds and related objects 
as organizing structures, and probabilistic methods, which study the evolution of probability 
densities. Geometric methods include the study of invariant manifolds, the theory of lobe 
dynamics 

QMM in two- (and some three-) dimensions, and the notions of finite-time hyper- 
bolic material lines Q and surfaces Js(] . The latter objects are often studied computationally 
via finite-time Lyapunov exponent (FTLE) fields (3, 9j. All of these geometric objects rep- 
resent transport barriers and in this way influence (mitigate) global transport. Probabilistic 
approaches include a study of almost-invariant sets and very recently, coherent sets 



[141 . Il5l |. For autonomous and time-dependent systems respectively, almost-invariant sets 
and coherent sets represent those regions in phase space which are minimally dispersed un- 
der the flow. Such regions provide an ordered skeleton often hidden in complicated flows. 
A recent comparison of the geometric and probabilistic approaches is given in [16] for the 
time-independent and time-periodic settings. 

The probabilistic methodologies provide important transport information that is often 
not well resolved by geometric techniques. Minimally dispersive regions need not be identi- 
fied by geometric approaches. For example, recent work [16[ has shown that regions enclosed 
by FTLE ridges need not represent maximal transport barriers. Several authors 



have noted other shortcomings of the FTLE-based approach: potential ambiguity in mul- 
tiple FTLE "ridges", ambiguity over flow duration for FTLE calculations, and a lack of 
correspondence between the strength of the ridge and the dispersal of mass across the ridge. 

Probabilistic techniques have also been shown to be valuable analysis tools for geophys- 
ical systems. In such systems, physical quantities are often used to determine transport 
barriers. For example, lines of constant sea surface height (as proxies for streamlines under 
the assumption of geostrophy) are commonly used to determine locations of rotational trap- 
ping regions such as anticyclonic eddies and gyres |19| , and maximum gradients of potential 
vorticity (PV) are used to determine "edges" of vortices in the stratosphere [2CH22|. In both 
of these geophysical settings, the use of physical quantities has been shown to be non-optimal 



in determining the location of transport barriers [2 



Probabilistic and transfer operator approaches |lCHl3l. 1 161 . 125M27I] have proven to be very 
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effective for autonomous systems. Initial progress has been made in the development of these 
techniques for time-dependent systems [141 . LL5I ] over infinite time horizons. In the present 
work, we focus on transport analysis of time-dependent systems over a finite period of time, 
significantly expanding on concepts introduced in [241 ] , a nd developing finite-time analogues 
of the time- asymptotic coherent set constructions in [H, [l5|. We develop methodologies to 
identify those regions of phase space that are minimally dispersive, or maximally coherent, 
under the flow, for a specific finite time interval. We demonstrate the efficacy of our approach 
on two examples: an idealized stratospheric flow and a flow obtained from assimilated data 
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sourced from the European Centre for Medium Range Weather forecasting. In the first 
example, we demonstrate that our new methodology easily detects an important dynamical 
separation of the domain; this separation is not clearly evident from an examination of the 
finite-time Lyapunov exponent (FTLE) field 0-0] • In the second example, we show that our 
new techniques can isolate the Antarctic polar vortex to high accuracy on a two-dimensional 



isentropic surface when compared to the commonly used potential vorticity criterion |20| . 
We also illustrate our technique directly in three dimensions, going beyond the capabilities 
of existing techniques to image the vortex in three dimensions. 

An outline of the paper is as follows. In section [2] we describe our setting and outline 
our main computational tool, the transfer operator (or Perron- Frobenius operator), and our 
numerical approximation approach. In section [3] we motivate and detail our new compu- 
tational approach. Section 13.21 describes the necessary computations and sections H] and [5] 
illustrate our new methodology via two case studies. 



2 Flows, coherent sets, and transfer operators 

Let M C R d be a compact smooth manifold and consider a time-dependent vector field 
f{z ) t) ) z e M, tel. Suppose that / is smooth enough for the existence of a flow map 
$(z,t,r) :MxKxR-) M, which describes the terminal location of an initial point z at 
time £, flowing for r time units. 

Given a base time t and a flow duration r, our motivation is to discover coherent pairs 
of subsets A t) A t+r C M such that r) « A t + T . More precisely, we will call A t) A t+r a 

(p , t-> t) — coherent pair if 

p„(A u A+r) := K A t n $(A+r, t + r; -r))/fi(A t ) > p , (1) 

and n(A t ) = /i(At +r ), where \i is a "reference" probability measure at time t. The measure 
H describes the mass distribution of the quantity we wish to study the transport of over the 
interval [t, t + r]; \i need not be invariant under the flow $. 

We are only interested in coherent pairs that remain coherent under small diffusive per- 
turbations of the flow; robust coherent pairs. Clearly, a (1, i, r)-coherent pair can be pro- 
duced by choosing an arbitrary A t and setting A t+T = $(A t ,£;r). However, such a pair 
may not be stable if some diffusion is added to the system. In a chaotic system, the set 
At+ T = &{A t) t]r) defined as above will experience stretching and folding, and for moderate 
to large r will become very thin and geometrically irregular. A small amount of diffusion 
will then easily eject many particles from A t+T) reducing the coherence ratio Pn(A t , A t + r ). 
The requirement that coherent pairs be robust under diffusive perturbations favors coherent 
sets that are geometrically regular; these robust, regular sets are more likely to be more 
dynamically meaningful than non-robust, irregular sets. 

Our basic tool for identifying sets satisfying (pQ) is the transfer (or Perron-Frobenius) 
operator V t , T > L x {M,i) O defined by 

V t , T f(z) := t + r; -r)) • | det D$(z, t + r; -r) | (2) 

where £ is normalized Lebesgue measure on M. If f(z) is a density of passive tracers at time 
t, Vt, T f(z) is the tracer density at time t + r induced by the flow In the autonomous 
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setting, almost-invariant sets were determined jlll . Il2l . |25| by thresholding eigenfunctions / 
of Vt, T (— ^Pt for all t) corresponding to positive eigenvalues A « 1: A = {f < c] or {/ > c}. 

The above calculations involved constructing a Perron- Frobenius operator for the action 
of $ on the entire domain M. In the time-dependent setting, we wish to study transport 
from X C M to a small neighborhood Y of £; r) C M. A global analysis would mean 
that X = Y = M and a transfer operator would be constructed for all of M. However, often 
one is interested in the situation where the domain of interest X is "open" and trajectories 
may leave X in a finite time (our numerical examples in Sections [4] and [5] illustrate this) . 
Moreover, the subset X may be very small in comparison to M. In such instances, there 
are great computational savings if the analysis can be carried out using a non-global Perron- 
Frobenius operator defined on X rather than M. Our new methodology allows precisely 
this and is a significant theoretical and numerical advance over existing transfer operator 
numerics. 

We now describe a numerical approximation of the action of Vt, T from a space of functions 
supported on X to a space of functions supported on Y. We subdivide the subsets X and Y 
into collections of sets {£>i, . . . , B m } and {Ci, • • • ? C n } respectively. We construct a finite- 
dimensional numerical approximation of the transfer operator 7\ r , using a modification of 
Ulam's method 



281 



where I is a normalized volume measure. Clearly, the the matrix P^(t) is row-stochastic by 
its construction. The value P( r \t)ij may be interpreted as the probability that a randomly 
chosen point in Bi has its image in Cj. We numerically estimate P^(t)ij by 

P (T) (t\j « #{r : Zi tr e Bi^iz^t-r) G C.-j/Q, (4) 

where z^ r , r = 1,...,Q are uniformly distributed test points in J3j(t) and $(z^ r ,t;r) is 
obtained via a numerical integration. 

The numerical discretization has the useful side-benefit of producing a discretization- 
induced diffusion with magnitude the order of the image of box diameters (see Lemma 2.2 
(29|). Ultimately, in Section [3T21 we will construct coherent sets by thresholding vectors in 
spjxsi j • • • j XB m } an d spjxci , • • • , Xc n }- This discretization limits the irregularity of possible 
coherent sets, and in practice, high regularity is observed. 



3 Coherent partitions 

For the remainder of the paper we set = P( r )(£)^j, fixing t and r. We set pi = fi(Bi), i = 
1, . . . , m and assume that pi > for all i = 1, . . . , n (if some sets Bi have zero reference 
measure, we remove them from our collection as there is no mass to be transported). Define 
q = pP to be the image probability vector on 7; we assume q > (if not, we remove 
sets Cj with qj = 0). The probability vector q defines a probability measure v on Y via 
v(Y f ) = Y^j=i Qj^O^^Cj) fo r measurable Y' C Y . We may think of the probability measure 
v as the discretized image of \i. 
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3.1 Problem setup 

To find the most coherent pair, we first try to partition X and Y as X\ U X 2 and Y\ U Y 2 in 
a particular way, where Xi,X 2 , Yi,F 2 all have measure approximately 1/2. This restriction 
will be relaxed later. Let Ji, 7 2 partition {1, . . . , m} and Ji, J 2 partition {1, . . . , n} and set 
X k = U ieIk Bi and F fc = U je j k Cj, k = 1,2. We desire: 

1. »(X k ) = Ziei k Pi « 1/2, Kn) = E iG j, % « 1/2, A; = 1, 2 

(the sets X\,X 2 and Yi,!^ partition X and F into two sets of roughly equal /i-mass 
and z/-mass respectively). 

2. p /i (X,,y,)^l,fc=l,2 

(this is a measure-theoretic way of saying that <&(X k , t\r) ~Y k , k = 1,2). 

3.2 Solution approach 

Introduce the inner products (x u x 2 ) p = Yi x i,i x 2,iPi and (y u y 2 ) q = Yj Vi^jQj- We form a 
normalized matrix = piPij/qj. The L resulting from this normalization of P ensures that 
1L — 1. We think of L as a transition matrix from the inner product space (]R m , -) p ) to the 
inner product space (]R n , (•, -) g ), which takes a uniform density on (]R m , -) p ) (representing 
the measure fi) to a uniform density on (]R n , -) q ) (representing the measure v). 

To describe 2-partitions of X and Y we consider vectors x G {±l} m , ?/ G {±l} n and define 
Xi = U, t= i^,X 2 = U, i= -i^^i = IWa,^ = U <w =-i^i- Thus > the partitions 
7i,/ 2 and Ji, J 2 are described by the parity of x and y respectively. We can write the 
condition ji{X k ) = Yi^i k Pi ~ 1 / 2 ^( Y k) = Qi ~ V 2 , fc = 1,2 as l) p |, < e 

for small e (the e is needed as it may be impossible to form finite collections of sets B { and 
Cj with measure exactly 1/2). 

Consider the problem: 

max{^L, y) q :xe {±l} m , y G {±l} n , \(x, l) p \, \(y, l) q \ < e} (5) 

for some small e. 
The objective 



{xL,y) q = 




(//(Xx n $(Y u t + r; -r)) + //(X 2 n $(F 2 , t + r; -r))) 
- (^(Xx n $(F 2 , * + r; -r)) + ^(X 2 n $(Y u t + r; -r))) 



= ^(X 1 )p fl (X 1 , Y l ) + n(X 2 ) Pfl (X 2 , Y 2 ) - ^(X 1 )p fl (X 1 , Y 2 ) - p(X 2 )p,(X 2 , Y t ). 

Thus, maximizing (xL, y) q is a very natural way to achieve our aim of finding partitions so 
that p^(X k ,Y k ) « = 1,2. The approximation in the above reasoning occurs because 
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Pii « p{B % n $(C J5 1 + r; _ T ))//i(Sifl 

The problem (j5J) is a difficult combinatorial problem; as a heuristic means of finding a 
good solution we relax the binary restriction on x and y and allow them to take on continuous 
values. We will interpret the values of x and y as "fuzzy inclusions" ; if Xi is very positive, 
then Bi is very likely to belong to Xl, and if Xi is very negative, then B { is very likely to 
belong to X 2 . Similarly for strong positivity or negativity of yi and inclusion of Bi in Y\ or 
Y 2 respectively. If the value of Xi or yi is near to zero, the fuzzy inclusion is less certain and 
we use an optimization in Algorithm ffl (Section 3.3) to determine where Bi belongs. 

As x and y can now float freely, we can set e = 0, and thus may insist that (x, l) p = 
(y, l) q — 0. When restricting x and y to be elements of {±l} m and {±l} n , we implicitly set 
the norms \\x\\ p = (x, x)]J 2 and \\y\\ q = (y, y) X J 2 to both be 1. Now that we let x and y freely 
float, we must include normalization terms in our objective. Thus, the relaxed problem is 

•&{pnit :<x ' 1> ' =(!/ ' 1> ' =0 } (6) 

We will use the optimal x and y to create our partition X\^X 2 and Y\,Y 2 via X\ = 

Ui: Xi >b B ii X 2 = U:^<A y i = Ui: yi >c C i> Y 2 = Oi:y i<c C h where b and c are chosen so 
that "}2 ieIk Pi ~ 1/2 X^ G j fc Qi ~ 1/2, fc = 1, 2. As an extension to our heuristic, we may also 
relax the condition that the measures of Xi ) X 2 and Yi,l2 are all approximately 1/2, and 
only enforce ji{Xk) = J2iei k Pi ~ Y2jej k Qi = v(Yk),k = 1,2. This would mean that while 
there is some flexibility in the choice of 6, the value c is a function of 6; see Algorithm [TJ 
We close this section with a lemma stating the solution to (|6]). 

Lemma 1. Let U p be an m x m diagonal matrix with p on the diagonal and U q be an n x n 
diagonal matrix with q on the diagonal Suppose that PP T is an irreducible matri^. The 
value of is o 2l the second largest singular value of U^ 2 PUq 1 ^ 2 , and the maximizing x 
and y in ^ are given by x = xllp 1 ^ 2 and y = yllq 1 ^ 2 , where x and y are the corresponding 
left and right singular vectors. 

Proof See appendix. □ 



3.3 Extraction of coherent pairs 

We now detail the procedure that extracts the coherent pairs from the vectors x and 

y identified in Lemma [U We create sets that are unions of boxes with x and y values above 
certain thresholds. Define X\{b) := Ur^>6^ anc ^ ^i( c ) := V]j- yj > c ^v ^ c ^ D e fi ne 

p(Xi(6),Yi(c)) = . V ;, ' ! -' :r ' Vi; " / " / ''. (7) 

Z-^i:BiCXi(b) Pi 

The quantity p measures the discretized coherence for the pair Xl(6), Yi(c). Our proce- 
dure to vary the thresholds b and c so as to select Xi(b) and Yi(c) with largest p value is 
summarized below: 

x If fi is absolutely continuous with a positive density that is Lipschitz on the interior of each then 
this error goes to zero with decreasing diameter of Bi and Cj] see Lemma 3.6 (30| . 
2 there exists a k such that (PP T ) k > 0. 
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Algorithm 1. 

1. Let 7/(6) = argmin c / G M|/i(Xi(6)) — z/(Yi(c'))|. This is to make v(Y\(d)) as close as 
possible to fi(Xi(b)). 

2. Set 6* = arg max p(X\ (6), Y]_(j](b))). The value ofb* is selected to maximize the coher- 
ence. 

3. Define A t = X 1 := X x {b*) and A t+T = Y 1 := Yi(r)(b*)). 

To obtain X 2 and F 2 , we define X 2 = X 2 (6*) := LL,<&* B i and Y 2 = Y 2 (r)(b*)) := 
U^<7j(&*) ^he complements of Xi and Yi in X and y, respectively. Thus, we select X\ 
and Yi to be the most coherent pair and define X 2 and Y 2 as their respective complements. 
One now should repeat Algorithm [1] with X 2 {b) := Ur^<6 ^ and ^( c ) := U^< c Cj, b, c G K. 
in place of Ai(6) and Yi(c) to search from "the negative end" of the vectors x and y ) possibly 
picking up a pair with higher coherence, and defining Xl, Yi as the complements of X 2 , Y 2 . 

4 Example 1: Idealized Stratospheric flow 

We consider the Hamiltonian system ^ = — ^ = ff where 

= c^y — UoLtanh(y/L) + A 3 UoL8ech 2 (y/L) cos(fcix) 
+ A 2 [/ ^sech 2 (?//L) cos(/c 2 x — a 2 i) + Ai[/ ^ sec h 2 (?//I/) cos(fciX — ait). 

This quasi-periodic system represents an idealized stratospheric flow in the northern or south- 
ern hemisphere. Rypina et al |3j| show that there is a time-varying jet core oscillating in a 
band around y — and three Rossby waves in each of the regions above and below the jet 
core. The parameters studied in [31] are chosen so that the jet core forms a complete trans- 
port barrier between the two Rossby wave regimes above and below it. We modify some of 
the parameters to remove the jet core band and allow transport between the two Rossby wave 
regimes. We expect that the two Rossby wave regimes will form time-dependent coherent 
sets because transport between the two regimes is considerably less than the transport within 
regimes. We set the parameters as follows: c 2 /Uo = 0.205, c 3 /[/ = 0.700, As = 0.2, A 2 = 0.4 
and Ai = 0.075, with the remaining parameters as stated in Rypina et al [31] . 

Our initial time is t = 20 days and our final time is t + r = 30 days. At our initial time 
we set X = S 1 x [—2.5, 2.5] Mm, where S 1 is a circle parameterised from to 6.3717T Mm, 
and subdivide X into a grid of m — 28200 identical boxes X = {Si, ... , B m }. This choice 
of m is sufficiently large to represent the dynamics to a good resolution. We compute an 
approximation of $(A, 20; 30) by uniformly distributing Q = 400 sample points in each grid 
box and numerically calculating $(^ r , 20; 30) using the standard Runge-Kutta method. The 
choice of Q is made so that over the flow duration, the image of boxes is well represented 
by the Q sample points per box. These Q x m image points are then covered by a grid of 
n = 34332 boxes {Ci, . . . , C n } of the same size as the i — 1, . . . , m. 

We set Y = \J^ =1 Cj : covering the approximate image of X. The transition matrix 

P = P 2 o°^ is computed using (|lj). 
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Figure 1: (a) The optimal vector x\ (b) the optimal vector y\ (c) Backward-time (blue) and 
Forward-time (red) FTLEs at t = 20 days computed with the flow time r = 10 days; (d) 
Backward-time (blue) and Forward-time (red) FTLEs at t = 30 days computed with the 
flow time r = 10 days. 



As the flow is area preserving, a natural reference measure /i is Lebesgue measure, which 
we normalize so that /jl(X) = 1. Thus, fi(Bi) — pi — i = 1, . . . , m and so (U p )u = 1/ra, 
i — 1, . . . , to. The vector q is constructed as q = pP. We compute the second largest singular 
value of Ftp 7 PU q 1 and the corresponding left and right singular vectors and thus determine 
x and y from Lemma [U The top two singular values were computed to be o\ = 1.0 and 
a 2 ~ 0.996. We expect x to determine coherent sets at time t = 20 days and y to determine 
coherent sets at time t + r = 30 days. Figure 1(a) and 1(b) illustrate the vectors x and y, 



which provide clear separations into red (positive) and green (mostly negative) regions. 

We apply the thresholding Algorithm CD to the vectors x and y to obtain the pairs (Xl, Yi) 
and (X 2 ,y 2 )i shown in Figures [2(a)] and |2(b)| To demonstrate that Yi w $(Xi,20; 10), we 



plot the latter set in Figure 2(c). When compared with Figure [2(b)] we see that there is very 
little leakage from Yi, just a few thin filaments. Similarly, Figures ]2(d)1and[2(b)1 compare Y 2 
and $(X 2 ,20; 10), again showing a small amount of leakage. This leakage is quantified by 
computing p(X u Y r ) w p(X 2 , Y 2 ) w 0.98. 

We compare our results with the attracting and repelling material lines computed via 
the finite-time Lyapunov exponent (FTLE) field [7] with the flow time r = 10. The ridges of 



the FTLE fields are commonly used to identify barriers to transport. Figures 1(c) and 1(d) 
present an overlay of forward- and backward-time FTLEs at t = 20 and t = 30, respectively. 
In this example, there are several FTLE ridges in the vicinity of the dominant transport 
barrier across the middle of the domain, and also several ridges far away from this barrier. 
The FTLE ridges do not crisply and unambiguously identify the dominant transport barrier 
shown in Figures [T(a)1 and [1(b) 



3 When determining Xi and Yi, Algorithm [T] produced values 6* « 0.0077 and r](b*) w 0.0005. 
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Figure 2: (a) The sets X\ (red) and X 2 (blue); (b) the sets Y\ (red) and Y 2 (blue); (c) the 
set 20; 10); (d) the set $(X 2 , 20; 10). 

5 Example 2: Stratospheric polar vortex as coherent 
sets 

In our second example, we use velocity fields obtained from the ECMWF Interim data set 
( |http: / / data.ecmwf .int / data/index.html ) . We focus on the stratosphere over the southern 
hemisphere south of 30 degrees latitude. In this region, there are strong persistent transport 
barriers to midlatitude mixing during the austral winter; these barriers give rise to the 
Antarctic polar vortex. We will apply our new methodology to the ECMWF vector fields in 
two and three dimensions to resolve the polar vortex as a coherent set. 

5.1 Two dimensions 

Our input data consists of two-dimensional velocity fields on a 121 x 240 element grid in the 
longitude and latitude directions, respectively. The ECMWF data provides updated velocity 
fields every 6 hours. The flow is initialised at September 1, 2008 on a 475K isentropic surface 
and we follow the flow until September 14. To a good approximation isentropic surfaces are 
close to invariant over a period about two weeks |32| . 

We set X = S 1 x [—90°, —30°], where S 1 is a circle parameterized from 0° to 360°. The 
domain X is initially subdivided into the grid boxes i = 1, . . . , to, where m = 13471 in 
this example. Based on the hydrostatic balance and the ideal gas law, we set the reference 
measure pi = Pr^ 7 a^ for alH = 1, . . . , to, where Pr^ is the pressure at the center point of Bi 
and &i is the area of box Bi. 

Using Q = 100 sample points z^ r , r = 1, . . . , Q uniformly distributed in each grid box 
Bi, i — 1, . . . , to we calculate an approximate image i; t)@ and cover this approximate 

4 We use the standard Runge-Kutta method with step size of 3/4 hours. Linear interpolation is used 
to evaluate the velocity vector of a tracer lying between the data grid points in the longitude-latitude 
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image with m = 14395 boxes {Ci, • • • 5 Cn} to produce the image domain Y. We construct 
P = Pr^ as described earlier using the same Q x m sample points. 

We compute x and y as described in Lemma [lj graphs of these vectors are shown in 
Figure [3] (Upper left and Upper right). Figure [3] (Lower left and Lower right) shows the 
result of Algorithm [lj extracting coherent sets A t and A t+r from the vectors x and y. We 
calculate the coherent ratio p^{A t) A t + T ) « 0.991, which means that 99.1% of the mass in A t 
(September 1, 2008) flows into A t + T (September 14, 2008), demonstrating a very high level 
of coherence. 



1 September 2008 



14 September 2008 




-0.01 -0.005 0.005 0.01 0.015 

1 September 2008 





Figure 3: [Upper left]: Graph of x (September 1, 2008). [Upper right]: Graph of y (September 14, 2008). 
[Lower left]: The red set represents the coherent set A t (September 1, 2008) obtained from Algorithm [TJ The 
green curve illustrates the vortex edge as estimated using P V. [Lower right] : As for [Lower left] at September 
14, 2008. 



To benchmark our new methodology, we will compare our result with a method commonly 
used in the atmospheric sciences to delimit the "edge" of the vortex. It has been recognized 
that during the winter a strong gradient of potential vorticity (PV) in the polar stratosphere 
is developed due to (1) strong mixing in the mid-latitudes (resulting from the breaking of 
Rossby waves emerging from the troposphere and breaking in the stratospheric "surf zone" 

coordinates. In the temporal direction the data is independently amnely interpolated. 
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j20| ) and (2) weak mixing in the vortex region. While potential vorticity depends only on the 
instantaneous vector field, potential vorticity is materially conserved for adiabatic, inviscid 
flow (both of which are good approximations in stratospheric flow over timescales of a week 
or two). Thus, PV may be viewed as a quantity derived from the Lagrangian specification 
of the flow and is therefore a meaningful comparator for these nonautonomous experiments 
who also use Lagrangian information. We used the method of Sobel et al. [221] to calculate 
a PV-based estimate of the vortex edge. The result is shown by the green curve in Figure [3] 
(Lower left and Lower right). Notating the area enclosed by the green curve at September 
1, 2008 by Af v and at September 14, 2008 by A^ T) we compute p^Af v ,A^ T ) « 0.984; 
98.4% of the mass in v flows into A^f T over the 13 day period. 

Our transfer operator methodology is clearly consistent with the accepted potential vor- 
ticity approach and in fact identifies a region that experiences slightly greater transport 
barriers across its boundary, indicated by the slightly larger coherence ratio: 99.1% versus 
98.4%. In the next section we apply our methodology in three dimensions to estimate the 
three-dimensional structure of the vortex. 

5.2 Three dimensions 

Strong transport barriers to midlatitude mixing in the southern hemisphere are also known 
to exist even in the full 3D case, where strong descent occurs near the edges of polar vortex at 
each pressure altitude j^sl, H|- I n principle, PV-based methods could be extended to three- 
dimensions by (i) slicing the three-dimensional region of interest into several nearby isentropic 
surfaces, (ii) applying the PV methodology on each individual isentropic surface to obtain 
an estimate of the vortex boundary on that surface, and (iii) stitching together these curves 
to form a reasonable two-dimensional surface, with the hope that the surface represents an 
estimate of the boundary of the three-dimensional vortex. This stitching together of several 
curves is a nontrivial computational task and complicated geometries may be missed by this 
relatively simple construction. The PV approach is likely to be more susceptible to noise 
than our direct approach because the computation of PV relies on estimates of derivatives 
of the velocity field (vorticity is the curl of the velocity field). Finally, such an approach 
would not utilise the full three-dimensional vector field, but rather a series of vector fields 
on isentropic surfaces. 

A key point of our new methodology is that it can easily applied in either two or three 
dimensions and works directly with the velocity fields to compute coherent regions with 
minimal external flux. 

We set X = S 1 x [-90°, -30°] x [50,70], where the third (vertical) component of this 
direct product is in units of hPa. The ECMWF data is again provided on a 240 x 121 grid in 
the longitude/latitude directions, and additionally at 7 pressure levels between 20 and 150 
hPa. We use the full 3D velocity field from the ECMWF reanalysis data. 

We subdivide X into a grid of m = 4116 x 8 = 32928 (longitude-latitude x pressure) 
boxes, where all boxes have the same area in the longitude-latitude directions and a "height" 
of (70 — 50) /8 = 20/8 hPa in the pressure direction. Following hydrostatic equilibrium 
considerations, we set the mass pi of box Bi to be proportional to the base area of Bi 
multiplied by the box "height" in hPa, and normalise so that Y^i % Vi = 1- We select 
Q = 250 sample points in each grid box, uniformly distributed in the longitude-latitude 
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direction and equally spaced in pressure direction. The Q xm images of these sample points 
are then covered by a grid of n = 51722 boxes. 

Repeating the approach of the two-dimensional study, the two largest singular values are 
computed to be G\ « 1.0 and a 2 ~ 0.9994. A slice along the uppermost pressure level (50 
hPa) of the optimal vectors x and y is shown in Figure @J 



50hPa, 1 September 2008 



SOhPa, 7 September 2008 




(a) 




(b) 



Figure 4: A slice of 3D vectors x and y along the 50hPa pressure level. 

Applying Algorithm [U we compute the coherent sets A t and A t + T shown in Figure [5] 
with p^Ai, A u ) « 0.9890. Figures 5 (a) | and 5(b)] show that at 1 September 2008, a compact 
central domain with nearly vertical sides is extracted by Algorithm [TJ Figure 5(e) shows 
that after 6 days of flow, this set is advected both upwards and downwards, and that this 
advection is not uniform over all latitudes. Figure [5(e)] and Figure 5(f)| (which gives a view 
from "below"), demonstrate that the upward flow occurs primarily near the centre of the 
vortex (high latitudes) , while the downward flow is concentrated around the periphery (lower 
latitudes) . A bowl-like shape is evident in Figure |5(f) showing a thin layer at the core of the 
coherent set at 7 September, descending toward the troposphere near the edge of coherent 
set. This observation agrees with the motion of ozone masses in the lower stratosphere, 
where the mass in the mixing zone around the mid-latitude slowly_moves downward and the 



mass in the vortex core moves within a thin stratospheric layer [33|, [34 



6 Conclusions 

We introduced a methodology for identifying minimally dispersive regions (coherent sets) 
in time-dependent flows over a finite period of time. Our approach directly used the time- 
dependent velocity fields to construct an ensemble description of the finite-time dynamics; 
the Perron- Frobenius (or transfer operator). The transport of mass is explicitly calculated 
in terms of a reference measure considered to be most appropriate for the application by the 
practitioner. Singular vector computations of matrix approximations of the Perron- Frobenius 
operator directly yielded images of the coherent sets; the left singular vector described the 
coherent region at the initial time and the right singular vector at the final time. Our 
methodology is the first systematic transfer operator approach for handling time-dependent 
systems over finite time durations. A particular feature of our approach is that one can focus 
on small subdomains of interest, rather than study the entire domain; this leads to major 
computational savings. 
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Figure 5: (a)-(c) show the optimal coherent sets at 1 September 2008 at different views. 
(d)-(e) show the optimal coherent sets at 7 September 2008 at different views. 



In our first case study we used this new technique to show that an idealized strato- 
spheric flow operates as two almost independent dynamical systems with a small amount 
of interaction across two Rossby wave regimes. Our second case study utilised reanalysed 
velocity data sourced from the European Centre for Medium Range Weather Forecasting 
(ECMWF) to estimate the location of the Southern polar vortex. Studying the dynamics 
on a two-dimensional isentropic surface, we found excellent agreement with traditional po- 
tential vorticity (PV) based approaches, and improved slightly over the PV methodology in 
terms of the coherence of the vortex. We also used the full three-dimensional velocity field 
to determine the vortex location in three dimensions, a computation not easily carried out 
with standard applications of the PV approach. 



A Proof of Lemma U 

We first show that the condition on y in ([6j) is unnecessary. 

max ( lf^4k ■ (x, l)v = o\ = max { ^ : (x, l) p = 1 = max : (x, l) p = 



(9) 
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with the maximizing y being y = xL. Setting y = xL, we see that (y, l) q = (xL, l) q = 
(x, 1L*) P = (x, l) p = (it is straightforward to check 1L* = 1; L* = P T ). Thus, since the 
maximizing y in ([9]) satisfies (y, 1) 9 = when (x, l) p = 0, we see that the value of ([6]) equals 
the value of the LHS of ([9]), and both ([6]) and ([9]) have the same maximizing x and y. 

We now convert the RHS of ([9]) to a maximization in the standard £ 2 norm by noting 
that (xi,x 2 ) p = (xiny 2 ,x 2 ny 2 ) 2 and (yuy 2 ) q = (yi^ 2 , y 2 ng /2 ) 2 . 

RHS of © = max [t^!k . (rf W 2 , in*/ 2 ), = ol = max . ( f p i/2) = \ , 

[ ||xny 2 || 2 p p J [ \\x\\ 2 J 

(10) 

1/2 

where we have made the substitution x = xllp . We claim that the leading singular value 
of lip 1 ^ 2 LIiy 2 = Yll^PUq 1 ^ 2 is 1, with corresponding left singular vector p 1 / 2 . 

To prove this claim, we show that 1 is the leading singular value of L with corresponding 
left singular vector 1 (where L is always considered as a linear mapping from (•, -) p to -) q ). 
Since 1L = 1 and 1L* = 1, one has ILL* = 1; also LL* is irreducible iff PP T is irreducible. 
By the Perron- Frobenius Theorem (eg. Thm 1.4 and 2.1 [35[), 1 is the largest real eigenvalue 
of LL*, and is simple; hence the largest singular value of L is 1 and the left and right singular 
vectors are 1 G M m and 1 GM n respectively. 

The result now follows from the Courant-Fischer theorem for symmetric matrices (see eg. 
Thm. 4.2.11 (36|), standard properties of singular vectors and the computation y = xL = 
xI\^ 2 L = (xllp 1 ^ 2 Lny 2 )IIg 1 ^ 2 = yllq 1 ^ 2 where y is the right singular vector of Hp 1 ^ 2 LH^ 2 
corresponding to a 2 . 
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